#include "cont_term.h"

int Cont_term::phase(int n) { return n % 2 ? -1 : 1; }

int Cont_term::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

int Cont_term::Tri_state(const State_Jacobi &state) {
  int unit = 1;
  int l12_2 = 2 * state.l12;
  int l3_2 = 2 * state.l3;
  if (Tri(l12_2, state.S12_2, state.J12_2) != 1) {
    return -1;
  }
  if (Tri(l3_2, unit, state.j3_2) != 1) {
    return -1;
  }
  if (Tri(state.J12_2, state.j3_2, state.J_2) != 1) {
    return -1;
  }
  if (Tri(state.T12_2, unit, state.T_2) != 1) {
    return -1;
  }
  return 1;
}

void Cont_term::init(int n_max) {
  this->pf_cont = 1.0 * LECs.LEC_CE * pow(hbarc, 6) /
                  (pow(F_pi, 4) * Cut_LambdaX * 2 * sqrt(3.0) * pow(Pi, 4));
  //cout << "good 1 ?" << endl;
  this->PFK_cont_cal(n_max);
  //cout << "good 11 ?" << endl;
}

double Cont_term::cal(const State_Jacobi &final_state,
                      const State_Jacobi &init_state) {
  // calculate contact term
  double value = 0.0;
  if ((Tri_state(final_state) != 1) || (Tri_state(init_state) != 1)) {
    cout << "can`t contain Tri condition in Cont_term.cal";
    exit(0);
    return 0.0;
  }
  // !!!! should n't have E conservation !!!!
  // if ((init_state.E != final_state.E))
  //  return 0.0;
  if ((init_state.J_2 != final_state.J_2) ||
      (init_state.T_2 != final_state.T_2))
    return 0.0;
  if ((init_state.l3 != 0) || (init_state.l12 != 0))
    return 0.0;
  if ((final_state.l3 != 0) || (final_state.l12 != 0))
    return 0.0;
  if ((final_state.S12_2 != init_state.S12_2))
    return 0.0;
  if ((init_state.S12_2 != init_state.J12_2) ||
      (final_state.S12_2 != final_state.J12_2))
    return 0.0;
  if ((init_state.j3_2 != 1) || (final_state.j3_2 != 1))
    return 0.0;

  if (init_state.T12_2 != final_state.T12_2)
    return 0.0;
  int phase_cont;
  phase_cont = phase(init_state.T12_2 / 2 - 1);

  /*int phase_test_t =
      init_state.n12 + init_state.n3 + final_state.n12 + final_state.n3;
  phase_cont = phase_cont * phase(phase_test_t);
  */

  double part1 = sixJ.quick(1, 1, init_state.T12_2, 1, 1, 2);
  // cout << "good ?" << endl;
  // cout << " n12 = " << init_state.n12 << "\t n3 = " << init_state.n3 << endl;
  // cout << " n12p = " << final_state.n12 << "\t n3p = " << final_state.n3
  //     << endl;

  double part2 = PFK_cont_vec[init_state.n12][init_state.n3] *
                 PFK_cont_vec[final_state.n12][final_state.n3];
  // cout << "part 2 = " << part2 << endl;
  value = this->pf_cont * phase_cont * part1 * part2;
  // cout<<"pre-fac = "<<pf_cont<<" phase_cont = "<<phase_cont<<"part1 =
  // "<<part1<<"part2 = "<<part2<<endl;
  return value;
}

void Cont_term::PFK_cont_cal(int n_max) {
  // calculate the PFK vector

  this->PFK_cont_vec.resize(n_max + 1);

  for (int n12 = 0; n12 <= n_max; n12++) {
    this->PFK_cont_vec[n12].resize(n_max + 1);
    for (int n3 = 0; n3 <= n_max; n3++) {
      // cout << "n12 = " << n12 << "\t n3 = " << n3 << endl;
      double part_2 = 0.0;
      for (int k12 = 0; k12 < wave.k_mesh; k12++) {
        double part_1 = 0.0;
        for (int k3 = 0; k3 < wave.k_mesh; k3++) {
          part_1 += wave.Pwave_vecwx[0][n3][0][k3] * wave.F_reg_vec[k3][k12];
          // part_1 +=
          // wave.Pwave_x[0][k3]*wave.Pwave_w[k3]*wave.Pwave_vec[n3][0][k3];//*wave.F_reg_vec[k3][k12];
        }
        // part_2 += part_1 * wave.Pwave_vecwx[0][n12][0][k12];
        part_2 += part_1 * wave.Pwave_x[0][k12] * wave.Pwave_w[k12] *
                  wave.Pwave_vec[n12][0][k12];
        part_1 = 0.0;
      }
      PFK_cont_vec[n12][n3] = part_2;
      part_2 = 0.0;
    }
  }
  // cout << "?? good ??" << endl;
}

void Cont_term::build() {

  int ch_size = jaco.channel.size();
  // cout<<"\t !! ch_size = "<<ch_size<<endl;
  this->vec.resize(ch_size);

  Pgs_Bar pgs;
  pgs.init(ch_size, 30);
  std::cout << "Begin to build contact term:" << std::endl;
  for (int i = 0; i < ch_size; i++) {
    pgs.PrintBar(i);
    int size = jaco.state[i].size();
    vec[i].resize(size);
    for (int x = 0; x < size; x++) {
      vec[i][x].resize(size);
      for (int y = 0; y < size; y++) {
        vec[i][x][y] = 0.0;
      }
      for (int y = 0; y <= x; y++) {
        int flag = -1;
        if (jaco.state[i][x].F > 0 && jaco.state[i][y].I > 0)
          flag = 1;
        if (jaco.state[i][y].F > 0 && jaco.state[i][x].I > 0)
          flag = 1;
        if (flag < 0)
          continue;
        double val = this->cal(jaco.state[i][x], jaco.state[i][y]);
        vec[i][x][y] = val;
        vec[i][y][x] = val;
        // cout << "channel : " << i << "\t x : " << x << "\t y : " << y
        //      << "\t val = " << val << endl;
      }
    }
  }
  pgs.EndBar();
  std::cout << "End to build contact term:" << std::endl;
}
void Cont_term::print(int ch_num) {
  int size = vec[ch_num].size();
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      if (vec[ch_num][i][j] != 0.0)
        cout << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
    }
  }
}
void Cont_term::print() {
  for (int ch_num = 0; ch_num < vec.size(); ch_num++) {
    int size = vec[ch_num].size();
    cout << "channel num : " << ch_num << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (vec[ch_num][i][j] != 0.0)
          cout << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
      }
    }
  }
}
void Cont_term::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  for (int ch_num = 0; ch_num < vec.size(); ch_num++) {
    int size = vec[ch_num].size();
    file << "channel num : " << ch_num << endl;
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {
        if (vec[ch_num][i][j] != 0.0)
          file << i << "\t" << j << "\t" << vec[ch_num][i][j] << endl;
      }
    }
  }
  file.close();
}
